Equation of state of an interacting Bose gas at finite temperature: 
a Path Integral Monte Carlo study 



O 
O 

(N 



(N 



5-1 

43 



C 

o 
o 



o 

O 

-i— > 



T3 

o 
o 



X 



S. PilatiW, K. SakW 2 ), J. Boronat( 2 \ J. Casulleras( 2 \ and S. GiorginiM 
W Dipartimento di Fisica, Universita di Trento and CRS-BEC INFM, 1-38050 Povo, Italy 
^ Departament de Fisica i Enginyeria Nuclear, Campus Nord B4-B5, 
Universitat Politecnica de Catalunya, E-08034 Barcelona, Spain 
(Dated: February 5, 2008) 

By using exact Path Integral Monte Carlo methods we calculate the equation of state of an inter- 
acting Bose gas as a function of temperature both below and above the superfiuid transition. The 
universal character of the equation of state for dilute systems and low temperatures is investigated 
by modeling the interatomic interactions using different repulsive potentials corresponding to the 
same s-wave scattering length. The results obtained for the energy and the pressure are compared 
to the virial expansion for temperatures larger than the critical temperature. At very low temper- 
atures we find agreement with the ground-state energy calculated using the diffusion Monte Carlo 
method. 

PACS numbers: 



I. INTRODUCTION 

In the last decade, after the first realization of Bose- 
Einstein condensation (BEC) in dilute systems of alkali 
atoms pj , the experimental and theoretical investigation 
of quantum degenerate gases has become one of the most 
active and fast developing fields in atomic, molecular and 
condensed matter physics The effect of interatomic 
interactions on the properties of ultracold Bose gases has 
been the subject of a deep and extensive research activ- 
ity. As for the theoretical side, the problem has been 
addressed from many different perspectives, both at zero 
and finite temperature, focusing on dynamical or equi- 
librium properties, in various dimensionalities and geo- 
metrical configurations both homogeneous and inhomo- 
geneous. Many different methods have also been used, 
from simple mean-field approaches to more advanced and 
essentially exact quantum Monte Carlo techniques 0, Q . 
In particular, the Path Integral Monte Carlo (PIMC) 
method allows one to calculate, for a given interatomic 
potential, the equilibrium properties of a bosonic system 
at finite temperature essentially without any approxima- 
tion. 

The PIMC technique has been applied in the context 
of ultracold Bose gases to investigate various thermody- 
namic properties in harmonically confined systems both 
in three |4| and lower dimensions [||, and for a detailed 
study of the critical behavior and of the superfiuid transi- 
tion temperature in three- [(J Q and two-dimensional Q 
homogeneous systems. 

In the present paper, we report on the results of 
a PIMC calculation of the equation of state of a 
three-dimensional homogeneous Bose gas as a func- 
tion of temperature and for different values of the in- 
teraction strength. The main focus is on the uni- 
versal character exhibited by the equation of state in 
terms of the reduced temperature T/Tj, , where Tj? = 
(2ttTl 2 /mk B )[n/((3/2)} 2 / 3 is the BEC transition temper- 
ature of an ideal gas of particles of mass m and number 



density n, and of the gas parameter na 3 , incorporating 
the effects of interatomic interactions at low tempera- 
tures through the s-wave scattering length a. We con- 
sider different repulsive model potentials (hard sphere, 
soft sphere and negative-power potential) and we explic- 
itly show the universal behavior of the energy per particle 
and pressure, both below and above the transition tem- 
perature, if the gas parameter is small enough. At low 
temperatures, we compare the calculated energy per par- 
ticle with the results of a Diffusion Monte Carlo (DMC) 
study carried out at T = 0|9j, and at high temperatures 
with the virial expansion of an interacting gas. We be- 
lieve that the present microscopic calculation could serve 
as a reference study for investigations of the thermody- 
namic properties of interacting Bose gases. 

The structure of the paper is as follows. In Sec. [D] 
we give a brief overview of the PIMC method of which 
we use two different implementations: a fourth-order ex- 
tension of the Trotter primitive approximation (for the 
negative-power potential) and the pair-product approx- 
imation (for the hard- and soft-sphere potential). In 
Sec. IIIII we discuss the results for the energy per par- 
ticle and gas pressure as a function of temperature and 
interaction strength. Finally in Sec. lIVI we draw our con- 
clusions. 



II. METHOD 

We consider a system of N particles described by the 
following Hamiltonian 



S 2 N 



^ = -^5> 2 + £nh-r,|) 



(1) 



i<j 



with different models for the spherical two-body inter- 
atomic potential V(r): 
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1) Hard-sphere (HS) potential, denned by 

V HS (r\ - / +°° ( r < fl ) > 

V ^ r >-\ (r > a) , 



(2) 



for which the hard-sphere diameter a corresponds to the 

s-wave scattering length. 

2) Soft-sphere (SS) potential, defined by 



V bb (r) 



V (r<i? ) 
(r>i? ) 



(3) 



with Vb > 0. In this case the scattering length is given 
by 



a = Rq 



tanh(Xo-Ro) 
Ko~R~o 



(4) 



with Kq = Vom/h 2 . For finite Vq one has always Rq > a, 
while for Vq — > +oo the SS potential coincides with the 
HS one with Rq = a. In the present calculation the range 
of the SS potential is kept fixed to the value Rq = 5a and 
the height Vq is determined to give the desired value of 
a. 

3) Negative-power (NP) potential, defined by 



V NP {r) =a/r p , 



(5) 



with a > and the integer p > 3. For this potential the 
scattering length is given by [ljj 



/ 2ma/Ti 2 



1/(p-2) 



r[(p-3)/(p-2)] 
r[(p-i)/(p-2)] ' 



(6) 



where r(x) is the Gamma function. In the present 
calculation we use p — 9, which yields a — 
(2ma/49^ 2 ) 1 / 7 r(6/7)/r(8/7). 

The universal regime in the plane na 3 -T/T® is ana- 
lyzed by performing PIMC simulations using the above 
three potentials with the same value for the gas param- 
eter na 3 . 

The partition function Z of a bosonic system with in- 
verse temperature (3 — (/cpT) -1 is defined as the trace 
over all states of the density matrix p = e~@ H properly 
symmetrized. The partition function satisfies the convo- 
lution equation 

z = / dRp(R,PR,/3) = / dR ( 7 ) 

x J cZR 2 ... J rfRjvfP(R-) R-2, r)...p(Rj 



PR, 



where r = /3/M, R collectively denotes the position vec- 
tors R = (ri, r2, rjv), PR denotes the position vectors 
with permuted labels PR = (fp(i) , r P(2) > r P(N)) an d 
the sum extends over the iV! permutations of JV parti- 
cles. In a PIMC calculation, one makes use of suitable 
approximations for the density matrix p(R, R', r) at the 
higher temperature 1/r in Eq. (JJJ and performs the 



multidimensional integration over R, R2,...,Rm as well 
as the sum over permutations P by Monte Carlo sam- 
pling 11]. The statistical expectation value of a given 
operator 0(H), 



p 



rfRO(R)p(R,PR,/3) , 



(8) 



is calculated by generating stochastically a set of config- 
urations {Ri}, sampled from a probability density pro- 
portional to the symmetrized density matrix, and then 
by averaging over the set of values {0(R,)}. 

Various approximations have been used for the density 
matrix at the high effective temperature M/ (3. In a first 
approach, one relies on the exact operator formula 



,-r(T+V) + ^-[T,V] =e _ T T e -7 



(9) 



and approximates it in the limit r — ► 0. The lowest order 
is known as primitive approximation (PA) 



-r(T+V) _ -tV/2-tT-tV/2 



+ 0(r 3 ) 



(10) 



and generate results for the energy with a quadratic de- 
pendence on t. The convergence to the exact result is 
guaranteed by the Trotter formula |12| , 



e~< t+ v)= lim (e- T 

M->oo V 



T e -rV 



(11) 



but, from a practical point of view, PA is not accurate 
enough for studying systems at temperatures below the 
superfluid transition |13|. In order to improve the accu- 
racy of this approach, we have calculated the properties 
of the gas with the NP potential by means of a higher- 
order scheme based on the operatorial decompositions 
proposed by Chin 14]. Chin's action (CA) is accurate 
to fourth-order in r but allows for a practical sixth-order 
dependence by adjusting some free parameters of the de- 
composition |15|. 

An alternative approximation for the high tempera- 
ture density matrix is based on the pair-product ansatz 

(ppa) na 

pCR,R^)=f]>(r < y < ,T)n p ^44 - w 



In the above equation p\ is the single-particle ideal-gas 
density matrix 



Pi(r l; r-,r) 



/ m 



\2irh 2 ', 



Y 2 e -( ri -r^ m /(2n 2 r) ) (13) 



and p re i is the two-body density matrix of the interacting 
system, which depends on the relative coordinates = 
r ■ — r' , divided by the corresponding 



Yi - Tj and r '.j 
ideal-gas term 



Prel ( r y i r ij ' 



r) 



3/2 



-y m /(4h 2 T) 



(14) 
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It can be shown [13 that PPA, Eq. CU, is more accurate 
than the simple PA, Eq. (|10|) . especially when the tem- 
perature of the system is very low and the interactions 
are of hard-core type. We have used the PPA approach 
for the simulations with the HS and SS potentials which, 
in fact, can not be strictly used in the first approach due 
to their discontinuous character. 

The two-body density matrix at the inverse tempera- 
ture r, p re ;(r, r', t), can be calculated for a given spheri- 
cal potential V(r) using the partial- wave decomposition 

p rei (r,r',r) = — £(2* + l)P £ (cos 9) (15) 

1=0 



where Pe(x) is the Legendre polynomial of order I and 9 
is the angle between r and r'. The functions Rk,e(r) are 
solutions of the radial Schrodinger equation 

h 2 



dr 2 



k.f 



2dR, 
r dr 

h 2 k 2 



£(£ + l) 



R 



m 



+ V(r)R k , e 
with the asymptotic behavior 



k,t 



2 sin(fcr-£7r/2 + <5 £ ) 
7r r 



(16) 



(17) 



holding for distances r » iij, where i?o is the range of 
the potential. The phase shift 5i of the partial wave of 
order i is determined from the solution of Eq. (|16(l for 
the given interatomic potential V(r). 

For the HS potential a simple analytical approxima- 
tion of the high-temperature two-body density matrix 
due to Cao and Berne 0] has been proven to be highly 
accurate. The Cao-Berne density matrix (r,r',r) is 
obtained using the large momentum expansion of the HS 
phase shift Si ~ ~ka + £n/2 and the large momentum ex- 
pansion of the solutions of the Schodinger equation l|16(l 
Rk,i(r) ~ \/2/7rsin[fc(r — a)]/r. This yields the result 



p%f(r,r',T) = 1 q(r + r') - a 2 

Prel( r > r '> T ) rr ' 

v -[rr'+a 2 -a(r+r')](l+co&9)m/{2h 2 T) 



(18) 



In the case of the SS potential, we calculate numer- 
ically prei(r, r', r) from Eqs. (fT^|l - (fT7jl . In the case of 
the HS potential, we use both the density matrix de- 
termined numerically and the Cao-Berne approximation 
[Eq. I|19fl ]. We have verified that in all cases the two 
procedures give indistinguishable results for the HS in- 
teraction within the present statistical uncertainty. 



III. RESULTS 

PIMC simulations have been carried out for a Bose gas 
with periodic boundary conditions and with TV ranging 



-0.1 

tr° -0.2 

m 

-0.3 

I- 

EQ 
it 

c -0.4 
-0.5 



CO 



LU 



-0.6 



-0.7 



I I I I I I 


1 1 1 1 






- /* 




/ 

■ i 






- 


6 / 




w 


na 3 =10" 6 ; 

i i i i 



10 



T/T„ 



FIG. 1: (color online). Energy per particle and pressure of a 
Bose gas in the normal phase as a function of temperature. 
The gas parameter is na 3 = 10 -6 . Solid symbols (blue online) 
refer to E/N - 3k B T/2: HS potential (circles), SS potential 
(diamonds). Open symbols (red online) refer to P/n — ksT: 
HS potential (circles), SS potential (diamonds). Statistical 
error bars are smaller than the size of the symbols. The virial 
expansion 1201 is represented by lines (blue online): HS po- 
tential (solid line), SS potential (long-dashed line). The virial 
expansion 1191 is represented by lines (red online) : HS poten- 
tial (short-dashed line), SS potential (dotted line). 



from 64 to 1024. In all the calculations finite size ef- 
fects have been checked to be smaller than the reported 
statistical uncertainty. 



A. Normal phase 



At high temperatures nX^ 1, where At = 

2nh 2 /niksT is the thermal wave length, the equation 
of state of the gas can be calculated from the virial ex- 
pansion 



PV 

Nk B T 



l + a 2 (T)nA T + 



(19) 



where we considered only the contribution arising from 
the second virial coefficient 02 (T). The corresponding 
virial expansion of the energy per particle can be calcu- 
lated using standard thermodynamic relations and one 
finds [13 



E 
N 



-k B T (1 + a 2 {T)n\^ + 



(20) 



For a non-interacting Bose gas the second virial coeffi- 
cient can be promptly calculated with the result aE] = 
— l/v2^, determined by statistical effects. For a gas of 
particles interacting through a repulsive interatomic po- 
tential 02 (T) can be calculated through a summation over 
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FIG. 2: (color online). Energy per particle and pressure of a 
Bose gas in the normal phase as a function of the tempera- 
ture for na 3 = 10 -4 . Same notation as in Fig. Q except for 
the results of the energy for the NP potential, shown here as 
squares (blue online). 
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FIG. 3: (color online). Energy per particle and pressure of a 
Bose gas in the normal phase as a function of the temperature 
for na s — 10 -2 . Same notation as in Figs. and [5] 



partial waves 0] 



a 2 (T) 



£=0,2,4,... 



dke -^/mk B T 



dS e (k) 
dk 



(21) 



For bosons, the sum in Eq. (|21|l only includes even par- 
tial waves. The £-th partial- wave phase shift 5i(k) in 
the above equation is obtained from the solution of the 
Schrodinger equation (|16[) for the given potential V(r) 



with the boundary condition (|17|l . If the thermal wave 
length is much larger than the range of the potential, 
At ^> -Ro, one obtains, to lowest order, the following 
result 



a 2 (T)-a° 2 = 2— 

AT 



(22) 



which only depends on the scattering length a. 

In Figs. 11131 we present the PIMC results for the energy 
per particle E/N and pressure P in the normal phase 
(T > T c ) for three different values of the gas parameter: 
na 3 = 10~ 6 (Fig.EJ), na 3 = 1(T 4 (Fig-EJ) and no 3 = 
(Fig. |2| • In order to emphasize the deviations from the 
classical results we plot the quantities E /N—3ksT/2 and 
P/n — ksT. We also plot the corresponding virial expan- 
sions from Eqs. (|19I) . I|2U|) with the second virial coefficient 
a2 (T) calculated using Eq. H21|) for the HS and SS poten- 
tials. For the smallest value of the interaction strength, 
na 3 = 10 -6 , we find very good agreement between the 
HS and SS results and with the virial expansions down 
to temperatures close to the transition temperature. At 
na 3 — 10~ 4 , the virial expansion still provides a good ap- 
proximation in the whole temperature regime and devia- 
tions are found only at the lowest temperatures T ~ T c °. 
On the other hand, universality is maintained only for 
low T since differences between the HS and SS potentials 
start to become visible for temperatures T/T® > 4. Fi- 
nally, for the largest interaction strength na 3 = 10 -2 , the 
universal behavior fixed by the scattering length a is lost 
in the whole temperature range. For the HS potential, 
agreement with the virial expansion is found only at the 
largest temperature. Concerning the results for the NP 
potential obtained with the CA approximation, notice 
that already for na 3 — 10~ 4 the statistical uncertainty is 
significantly larger than the one corresponding to results 
for the HS and SS potentials obtained using the PPA. 
This fact is due to the large separation in scale between 
the range of interactions and the mean interparticle dis- 
tance which occurs in dilute systems. For very small 
values of the gas parameter na 3 the algorithm based on 
the PPA, which is constructed from the exact solution of 
the two-body problem, converges much faster than the 
one based on the CA. For example, at T = 2T® and 
na 3 — 10~ 4 , the calculation for the NP potential has 
been performed using up to M = 200 beads in contrast 
to only M = 12 in the case of the HS and SS poten- 
tials. For the smallest value of the interaction strength, 
na 3 = 10 -6 , and especially for temperatures below the 
transition temperature, the calculation using the CA ap- 
proach becomes much more computationally demanding 
due to the large number of beads required. 



B. Superfluid phase 

The determination of the transition temperature from 
the equation of state is very delicate as it involves a slight 
change of the energy vs. T dependence at T c . Its precise 
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FIG. 4: (color online). Energy per particle and pressure of 
a Bose gas in the superfluid phase as a function of temper- 
ature. The gas parameter is na 3 — 10 -6 . Solid symbols 
(blue online) refer to E/N: HS potential (circles), SS poten- 
tial (diamonds). Open symbols (red online) refer to P/n: HS 
potential (circles), SS potential (diamonds). Statistical error 
bars are smaller than the size of the symbols. The horizontal 
bar (green online) corresponds to the ground-state energy per 
particle Eo/N of a HS gas calculated using DMC. The lines 
correspond to a non-interacting gas: the solid line (blue on- 
line) refers to the energy per particle and the dashed line (red 
online) to the pressure. 



FIG. 6: (color online). Energy per particle and pressure of 
a Bose gas in the superfluid phase as a function of the tem- 
perature for na 3 = 10 -2 . Same notation as in Figs. 2] and 
□ 



determination would require the calculation of the spe- 
cific heat. Other observables, like the superfluid density 
and the condensate fraction, give a direct signature of 
the transition The most reliable results so far give 

the following shift of the transition temperature T c 



(1.29±0.05)(na d ) 



3\l/3 



(23) 
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FIG. 5: (color online). Energy per particle and pressure of a 
Bose gas in the superfluid phase as a function of the temper- 
ature for na 3 = 10 -4 . Same notation as in Fig. [IJ except for 
the results of the energy for the NP potential, shown here as 
squares (blue online). 



holding for very small values of the gas parameter. For 
the values of na 3 used in the present study, the above 
equation yields estimates of T c ranging from 1% to 30% 
above T°. We are not interested here in the calculation 
of T c and the focus is only on the precise determination 
of the equation of state. 

In Figs. 14161 we show results for E/N and P/n in the 
superfluid phase (T < T c ) for the three values of the gas 
parameter: na 3 = 1CT 6 (Fig. EJ, na 3 = 1CT 4 (Fig. EJ) 
and na 3 = 10~ 2 (Fig. 0). To emphasize the effects of 
interactions in Figs. 0] and [S] we also plot the energy 
per particle and the pressure of an ideal Bose gas. For 
na 3 = 1CP 6 and na 3 = 10~ 4 we see very small differ- 
ences between the results of the HS and SS potentials 
(notice the large statistical uncertainty in the results of 
the NP potential at na 3 = 10~ 4 ). For the largest value 
of na 3 we still find good agreement between the results 
of the HS and NP potential, while significant differences 
are found between the HS and SS potential. At very low 
temperatures, the PIMC results agree with the ground- 
state energy per particle Eq/N obtained for the HS po- 
tential using the DMC method 0. For the three val- 
ues of the gas parameter used in the present study the 
results are (in units of fc B T c °): E /N = 0.01905(2) at 
na 3 = I0~ 6 , E /N = 0.09185(8) at na 3 = 10~ 4 , and 
Eq/N = 0.5840(4) at no 3 = 10" 2 . 
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The PIMC results for E/N and P/n obtained using 
the HS and SS potential at the various temperatures 
and for the three values of na 3 are reported in Tables [i] 
ITT1 Finite size effects are relevant only in the simula- 
tions performed at T = because of the vicinity to 
the critical point. For example, in the case of the HS 
potential at na 3 = 10~ 4 we find for this temperature 
the result: E/{Nk B T°) = 0.935(7) with N = 64 and 
E/(Nk B T°) = 0.891(13) with N = 1024. 



the gas parameter is pointed out by performing simula- 
tions with different interatomic model potentials. Above 
the transition temperature we compare our results for en- 
ergy and pressure with the high-temperature expansion 
based on the second virial coefficient. The inclusion of 
tables for both energy and pressure is intended as a refer- 
ence for future studies of the thermodynamic properties 
of interacting Bose gases. 



IV. CONCLUSIONS 

In conclusion, we have carried out using exact PIMC 
methods a precision calculation of the equation of state 
of an interacting Bose g CIS CIS cl function of temperature 
and interaction strength. The universal character of the 
equation of state at low temperatures and small values of 
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TABLE I: Energy per particle E/N (in units of fcsT c °) for the HS and SS potential and different values of the gas parameter 
na 3 = 10~ 6 , 10~ 4 , 10~ 2 . In parenthesis we give the statistical errors. 
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0.856(13) 


0.885(13) 


1.00 


0.816(10) 


0.819(6) 


0.891(13) 


0.902(15) 


1.398(10) 


1.431(18) 


2.00 


2.540(6) 


2.541(5) 


2.665(13) 


2.675(8) 


3.328(4) 


3.090(3) 


4.00 


5.694(5) 


5.692(5) 


5.818(7) 


5.821(5) 


6.458(3) 


6.196(3) 


10.0 


14.817(5) 


14.821(5) 


14.956(4) 


14.943(3) 


15.527(6) 


15.312(4) 


20.0 


29.885(5) 


29.880(5) 


30.023(4) 


29.989(2) 


30.615(11) 


30.374(4) 



TABLE II: Pressure P/n (in units of UbT°) for the HS and SS potential and different values of the gas parameter na 3 = 10 
10~ 4 , 10~ 2 . In parenthesis we give the statistical errors. 









P/n 


(k B T°) 








na 3 = 10~ 


6 


na? 


= 10- 4 


na 3 


= 10~ 2 




HS 


SS 


HS 


SS 


HS 


SS 


0.10 


0.021(1) 


0.021(1) 


0.095(1) 


0.094(1) 


0.676(6) 


0.514(1) 


0.25 


0.037(2) 


0.035(1) 


0.107(2) 


0.106(1) 


0.680(4) 


0.520(2) 


0.50 


0.116(3) 


0.114(2) 


0.183(4) 


0.180(4) 


0.724(4) 


0.582(3) 


0.75 


0.285(5) 


0.283(5) 


0.362(9) 


0.353(9) 


0.890(8) 


0.761(8) 


1.00 


0.556(7) 


0.558(4) 


0.648(10) 


0.651(10) 


1.276(8) 


1.113(12) 


2.00 


1.708(4) 


1.706(3) 


1.835(9) 


1.832(5) 


2.572(3) 


2.223(3) 


4.00 


3.808(4) 


3.807(3) 


3.937(5) 


3.924(4) 


4.716(4) 


4.308(3) 


10.0 


9.891(4) 


9.893(3) 


10.028(3) 


9.995(2) 


10.954(8) 


10.378(3) 


20.0 


19.936(3) 


19.932(3) 


20.075(4) 


20.024(2) 


21.335(10) 


20.416(3) 



